16S Phyloseq

Load required packages.

  require(phylosmith)

Loading required package: phylosmith

Load necessary files from Dada2 processing script.

load("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/16S/DADA2outputs/dada2_outputs_16S_BeeBread_wtree.RData")

Read in metadata

metadata<-read.csv("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/BB_sample-metadata.txt", sep = "\t")
rownames(metadata)<-metadata$X
metadata$Region_Site<-paste( metadata$Region, metadata$Site, sep = "")
metadata$Treatment_Site<-paste( metadata$Treatment, metadata$Site, sep = "")

Convert files to correct format for phyloseq object and then rename samples in OTU table to match the format of the metadata file. Create phyloseq object and root phylogenetic tree for usage with distance metrics (e.g., Unifrac).

#assign object type to merge into phyloseq object
md<-sample_data(metadata)
otu<-otu_table(seqtab.nochim, taxa_are_rows = F)
tax_tab<-tax_table(taxa)
#check sample names of both otu table and metadata
#remove .fastq.gz from otu table names
sample_names(otu)<-str_split(sample_names(otu), pattern = ".fastq.gz", simplify = T) [,1]
#since they match, you can create the phyloseq object
bb16S_orig<-phyloseq(md, otu, tax_tab, fitGTR$tree)

#root tree for distance metrics like unifrac
set.seed(11)
phy_tree(bb16S_orig)<-root(phy_tree(bb16S_orig), sample(taxa_names(bb16S_orig), 1), resolve.root = TRUE)
is.rooted(phy_tree(bb16S_orig))

Pre-processing. Rename OTU IDS to something more manageable. We save the ASCV sequences in case we need them later (e.g., BLASTN searches). Create Bacterial and Archaeal phyloseq objects, but first we remove Chloroplasts.

#Extract original IDs for future reference. 
seq_df_w_OTUID<-data.frame(OTUID = 1:ntaxa(bb16S_orig), Sequence = taxa_names(bb16S_orig))
taxa_names(bb16S_orig)<-paste("OTU_" , 1:ntaxa(bb16S_orig), sep = "")

#remove chloroplast
#60 chloroplast taxa
bb16S_orig_nc<-subset_taxa(bb16S_orig, Class != "Chloroplast")

#2257 bacterial taxa - 9 non-bacteria removed. 
bb16S_bac<-subset_taxa(bb16S_orig_nc, Kingdom == "Bacteria")
#9 archaea
bb16S_arch<-subset_taxa(bb16S_orig_nc, Kingdom == "Archaea")

Remove control samples as there is only as single replicate of each.

ps_wo_control<-subset_samples(bb16S_bac, Treatment != "none")

Examine rarefaction curves. Plateau in most samples under the 5k read mark.

rarecurve(otu_table(ps_wo_control), step=50, cex=0.5)

First, check sample depths for rarefaction. Minimum sample depth is a bit over 21k reads, so we will rarefy at 20k.

set.seed(11)
sort(sample_sums(ps_wo_control))
bb16S_bac_rarefy<-rarefy_even_depth(ps_wo_control, sample.size = 20000, rngseed = 14, trimOTUs = T)
#300 otus were removed during rarefaction. 

#To move forward, we can hellinger transform and go from there.
bb16S_bac_hellinger<-transform_sample_counts(ps_wo_control, function(x) sqrt(x / sum(x)))

#maximum liklihood point estimates
bb16S_bac_ML<- transform_sample_counts(ps_wo_control, function(x) x / sum(x))

#Exploratory analysis of Alpha diversity

Plots of Alpha Diversity by region, treatment, timepoint, and the combinations.

plot_richness(bb16S_bac_rarefy, x = "Region", measures = c("Shannon", "Observed", "Chao1"), color = "Region") + geom_boxplot() + ggtitle("Region")

plot_richness(bb16S_bac_rarefy, x = "Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Treatment")

plot_richness(bb16S_bac_rarefy, x = "Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Site") + geom_boxplot() + ggtitle("Site")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site by timepoint")

plot_richness(bb16S_bac_rarefy, x = "Region_Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Treatment") + geom_boxplot() + ggtitle("Region_Treatment")

plot_richness(bb16S_bac_rarefy, x = "Region_Site", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Site") + geom_boxplot() + ggtitle("Region_site")

plot_richness(bb16S_bac_rarefy, x = "Region_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Timepoint") + geom_boxplot() + ggtitle("Region_Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Treatment_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment_Timepoint") + geom_boxplot() + ggtitle("Treatment_Timepoint")

Based upon the site and region graphs, it appears that site contains more differences than region. There are no obvious differences among treatments. The site by timepoint graph shows what appears to be an interaction of before and after treatment and site. Some sites display lower diversity following treatment and others higher, with consistent trends across all sites though.

Statistical testing of alpha diversity metrics. Using ANOVA framework and checking for assumption of normality and correlation in predictors. region and site appear to be associated, so I will drop Region for this exploratory analysis. We can revisit if we decide Region is a more accurate predictor, but it appears that from the graphs above site provides more resolution than region. https://www.pluralsight.com/guides/testing-for-relationships-between-categorical-variables-using-the-chi-square-test

#create alpha diversity table and prep data to include grouping columns
richness_BB16S_bac<-estimate_richness(bb16S_bac_rarefy, measures = c("Shannon", "Observed", "Chao1"))
#add metadata
richness_BB16S_bac$Region<-sample_data(bb16S_bac_rarefy)$Region
richness_BB16S_bac$Site<-sample_data(bb16S_bac_rarefy)$Site
richness_BB16S_bac$Timepoint<-sample_data(bb16S_bac_rarefy)$Timepoint
richness_BB16S_bac$Treatment<-sample_data(bb16S_bac_rarefy)$Treatment
#add combo columns
richness_BB16S_bac$Region_Treatment<-sample_data(bb16S_bac_rarefy)$Region_Treatment
richness_BB16S_bac$Region_Site<-sample_data(bb16S_bac_rarefy)$Region_Site
richness_BB16S_bac$Region_Timepoint<-sample_data(bb16S_bac_rarefy)$Region_Timepoint
richness_BB16S_bac$Treatment_Timepoint<-sample_data(bb16S_bac_rarefy)$Treatment_Timepoint

#test of association between 
chisq.test(richness_BB16S_bac$Region, richness_BB16S_bac$Site, simulate.p.value = TRUE)

Tests of Shannon diversity - significant differences across sampling sites and interaction term of timepoint and site.

summary(mod)

                         Df Sum Sq Mean Sq F value   Pr(>F)    
Site                      5  6.008  1.2015  14.092 1.19e-07 ***
Treatment                 2  0.144  0.0722   0.847    0.437    
Timepoint                 1  0.130  0.1302   1.527    0.225    
Site:Treatment           10  1.383  0.1383   1.622    0.140    
Site:Timepoint            5  3.132  0.6265   7.347 7.82e-05 ***
Treatment:Timepoint       2  0.002  0.0011   0.013    0.987    
Site:Treatment:Timepoint 10  1.487  0.1487   1.743    0.108    
Residuals                36  3.070  0.0853                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Tests of Chao 1 - emphasizes importance of rare taxa. Site is significant as are the interactions with site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)
summary(mod)
#Site and the interactions with site are significant. 

TukeyHSD(mod, "Site")
plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site x Timepoint")
plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Site x Treatment")

Tests of Observed - unique taxa in a sample. Site is significant as are the interactions with site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Observed ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)
summary(mod)
#Summary of model shows site to be a significant predictor as well as the site interaction terms.
TukeyHSD(mod, "Site")
plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Observed"), color = "Timepoint") + geom_boxplot() + ggtitle("Site")
plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Observed"), color = "Treatment") + geom_boxplot() + ggtitle("Site")

Summary of alpha diversity: Rarefying at 20k reads sufficiently captures the diversity of each sample, as demonstrated by the plateau in our rarefaction curves. Site is the most useful predictor for determining differences in alpha diversity. Several interactions with site are present for all metrics of alpha diversity. No clear trends in alpha diversity differences. Some sites display shifts responding to treatment, while others do not. Some sites show increased alpha diversity at timepoint 2 and others the opposite.

#Multivaraite exploratory analysis

Info on strata argument. Briefly, this allows us to restrict permutations to within a group and control for an effect in an analogous way to a random effect.

https://stats.stackexchange.com/questions/350462/can-you-perform-a-permanova-analysis-on-nested-data/350504#350504?newreg=8e389dc237c54d87ab240486b8b0fe33 https://stats.stackexchange.com/questions/188519/adonis-in-vegan-order-of-variables-or-use-of-strata?noredirect=1&lq=1 https://stats.stackexchange.com/questions/459407/permanova-outputs-with-or-without-random-factor https://ichthyology.usm.edu/courses/multivariate/feb_7.pdf

Both with Bray-curtis and Unifrac distances.

#extract data from phyloseq object
sd_adonis<-data.frame(sample_data(bb16S_bac_ML))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(bb16S_bac_ML, method =  "bray")

#full model with no strata arguemnt. Essentiall the effect of each main effect without controlling for any potential groups. 
adonis(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations to time points. 
adonis(dist ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Timepoint, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis(dist ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

#another way to specify the same model using adonis2. Both provide same results as would be expected. No need to run again, merely a sanity check.
#perm<-how(nperm = 10000)
#setBlocks(perm)<-with(sd_adonis, Timepoint)
#adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = perm, by = "terms")

#weighted unifrac distance
dist.uf<-phyloseq::distance(bb16S_bac_ML, method =  "wunifrac")
#full model with no strata argument, essentially the effect of each main effect without controlling for any potential groups. 
adonis(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations to time points. 
adonis(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Timepoint, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

Based on results, site is the main determinant of microbiome composition (r2 of 0.51 and r2 of 0.38, BC and weighted unifrac, respectively). Time point is also significant but explains much less of the variation (r2 of 0.063 and r2 of 0.021, BC and weighted unifrac, respectively). Treatment is only significant for BC (r2 of 0.024). Both remain significant for BC dissimilarity when site and timepoint are controlled for using the strata argument. However, when site is controled for using weighted unifrac dissimilarity, the only significant predictor is the site x treatment interaction. Next, I will split by time point and examine the effect of each treatment while controlling for site effects using the strata argument (strata = site).

Take away: Site is the strongest predictor of community dissimilarity. This is true for both Bray-Curtis and weighted unifrac distances.

Split by time point as these are before and after treatment.

#split then extract info from each
#pre-treatment samples
bb16S_bac_ML_A<-subset_samples(bb16S_bac_ML, Timepoint == "A")

#extract data from phyloseq object for pre-treatment
sd_adonis<-data.frame(sample_data(bb16S_bac_ML_A))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(bb16S_bac_ML_A, method =  "bray")

#full models with no strata 
adonis(dist ~  Site * Treatment, data= sd_adonis, permutations = 10000)
adonis(dist ~   Treatment * Site, data= sd_adonis, permutations = 10000)
#models with strata to control for site differences
adonis(dist ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)
adonis(dist ~  Treatment * Site, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

#weighted unifrac distance
dist.uf<-phyloseq::distance(bb16S_bac_ML_A, method =  "wunifrac")
#full model with no strata argument, essentially the effect of each main effect without controlling for any potential groups. 
adonis(dist.uf ~  Site * Treatment,  data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis(dist.uf ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

Time point 1 (pre-treatment): Significant effect of Site (r2 of 0.73). Treatment is also a significant predictor (r2 of 0.04). This is something to note because this is supposed to be before treatment was applied. When permutations are restricted to within site, treatment and site remain significant. For weighted unifrac, treatment is insignificant (r2 of 0.025) whether site is controlled for or not.

#post-treatment samples
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")

#extract data from phyloseq object for post-treatment
sd_adonis<-data.frame(sample_data(bb16S_bac_ML_B))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(bb16S_bac_ML_B, method =  "bray")

#full models with no strata 
adonis(dist ~  Site * Treatment, data= sd_adonis, permutations = 10000)
adonis(dist ~   Treatment * Site, data= sd_adonis, permutations = 10000)
#models with strata to control for site differences
adonis(dist ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)
adonis(dist ~  Treatment * Site, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

#weighted unifrac distance
dist.uf<-phyloseq::distance(bb16S_bac_ML_A, method =  "wunifrac")
#full model with no strata argument, essentially the effect of each main effect without controlling for any potential groups. 
adonis(dist.uf ~  Site * Treatment,  data= sd_adonis, permutations = 10000)
#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis(dist.uf ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

Time point 2 (post-treatment): Site was again a significant predictor (r2 of 0.60) and treatment was borderline significant (r2 of 0.03). However, when the effect of site was controlled for by using strata = site, treatment becomes significant (r2 of 0.039). So when controlling for site differences, the treatments contribute to differences in microbiome structure. This was not the case for weighted unifrac dissimilarity. Treatment was not a significant predictor of community dissimilarity regardless of whether or not site differences were controlled for.

Take away: Site remains the most useful predictor of community dissimilarity. Treatment effects are only observed for non-phylogentic metrics of dissimilarity when the effect of site is controlled for.

#Visualize ordinations using two methods. CAP and NMDS

#significant global model so lets look at individual parameters
drop1(cca_mod, test = "perm", permutations = 100)

          Df    AIC      F   Pr(>F)   
<none>       109.21                   
Site       5 122.00 4.6923 0.009901 **
Timepoint  1 111.59 3.9549 0.009901 **
Treatment  2 107.50 1.0161 0.435644   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Constrained ordination using CAP

#Bray-Curtis
ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "bray", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. 

ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "wunifrac", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic()  + ggtitle("Weighted Unifrac CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. There are several outlier samples that ordinate far to the right on the CAP1 axis. 

Unconstrained ordination using NMDS

#Bray-Curtis
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "bray", k = 3, trymax = 1000)

Run 0 stress 0.09985802 
Run 1 stress 0.09985806 
... Procrustes: rmse 0.0001233703  max resid 0.0006955126 
... Similar to previous best
Run 2 stress 0.09993189 
... Procrustes: rmse 0.0032591  max resid 0.01894355 
Run 3 stress 0.1020888 
Run 4 stress 0.1038354 
Run 5 stress 0.1094871 
Run 6 stress 0.09974117 
... New best solution
... Procrustes: rmse 0.01677408  max resid 0.09915296 
Run 7 stress 0.1057786 
Run 8 stress 0.1053099 
Run 9 stress 0.1082896 
Run 10 stress 0.1096804 
Run 11 stress 0.1103444 
Run 12 stress 0.09985781 
... Procrustes: rmse 0.01663359  max resid 0.09817835 
Run 13 stress 0.0998581 
... Procrustes: rmse 0.01678077  max resid 0.0991411 
Run 14 stress 0.09974171 
... Procrustes: rmse 0.0002949829  max resid 0.001459841 
... Similar to previous best
Run 15 stress 0.1053098 
Run 16 stress 0.1030694 
Run 17 stress 0.1020602 
Run 18 stress 0.1040017 
Run 19 stress 0.103077 
Run 20 stress 0.09974164 
... Procrustes: rmse 0.0002844151  max resid 0.001386766 
... Similar to previous best
*** Solution reached

ord_nmds$stress

[1] 0.09974117

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. 

#weighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "wunifrac")

Run 0 stress 0.07903344 
Run 1 stress 0.07903323 
... New best solution
... Procrustes: rmse 8.386795e-05  max resid 0.0004041708 
... Similar to previous best
Run 2 stress 0.07903364 
... Procrustes: rmse 0.0003671032  max resid 0.002573766 
... Similar to previous best
Run 3 stress 0.07903332 
... Procrustes: rmse 3.997463e-05  max resid 0.0002766345 
... Similar to previous best
Run 4 stress 0.07903317 
... New best solution
... Procrustes: rmse 4.432156e-05  max resid 0.0002753805 
... Similar to previous best
Run 5 stress 0.07903311 
... New best solution
... Procrustes: rmse 8.683721e-05  max resid 0.0005822905 
... Similar to previous best
Run 6 stress 0.0790332 
... Procrustes: rmse 0.0001033843  max resid 0.0007320929 
... Similar to previous best
Run 7 stress 0.0790336 
... Procrustes: rmse 0.000240026  max resid 0.001703741 
... Similar to previous best
Run 8 stress 0.09216884 
Run 9 stress 0.07903312 
... Procrustes: rmse 9.946586e-06  max resid 4.859879e-05 
... Similar to previous best
Run 10 stress 0.07903333 
... Procrustes: rmse 0.0001615213  max resid 0.001119523 
... Similar to previous best
Run 11 stress 0.0790334 
... Procrustes: rmse 0.0001782196  max resid 0.001222783 
... Similar to previous best
Run 12 stress 0.07903312 
... Procrustes: rmse 3.152764e-05  max resid 0.0002108642 
... Similar to previous best
Run 13 stress 0.07903345 
... Procrustes: rmse 0.0001995874  max resid 0.001359699 
... Similar to previous best
Run 14 stress 0.07903358 
... Procrustes: rmse 0.0002079531  max resid 0.001281117 
... Similar to previous best
Run 15 stress 0.07903337 
... Procrustes: rmse 0.0001601027  max resid 0.001025945 
... Similar to previous best
Run 16 stress 0.07903343 
... Procrustes: rmse 0.0001961275  max resid 0.00133597 
... Similar to previous best
Run 17 stress 0.07903338 
... Procrustes: rmse 0.0001784725  max resid 0.00113875 
... Similar to previous best
Run 18 stress 0.07903363 
... Procrustes: rmse 0.0002480767  max resid 0.001720673 
... Similar to previous best
Run 19 stress 0.07903322 
... Procrustes: rmse 0.0001162361  max resid 0.0004919752 
... Similar to previous best
Run 20 stress 0.07903344 
... Procrustes: rmse 0.0001989235  max resid 0.001414217 
... Similar to previous best
*** Solution reached

ord_nmds$stress

[1] 0.07903311

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Weighted Unifrac NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#Indicator taxa using Baruta Algorith Specify function for Baruta that accepts a phyloseq object.

baruta_phyloseq<-function(physeq = "Phyloseq object", GroupingVar = "Grouping variable"){
  #prep data
  df<-data.frame((otu_table(physeq)))
  md<-data.frame(sample_data(physeq))
  tax<-data.frame(tax_table(physeq))
  rownames(tax)<-str_replace(rownames(tax), pattern = "[[.]]", replacement = "=")
  #run simper
 
                  boruta_layer<-Boruta(df , factor(md[,GroupingVar]), 
                      doTrace = 1, ntree = 1000, maxRuns = 1000) # adjust factor!!!!

                  B1<-data.frame(boruta_layer$finalDecision)  
                  B2<-data.frame(apply(boruta_layer$ImpHistory, 2, mean))
                  names(B1)<-"Layer_decision"
                  names(B2)<-"Layer_importance"
                  BB1<-merge(B1, B2, by="row.names", all.x=TRUE)
                  #sorting the results
                  BB1_sort<-BB1[order(BB1$Layer_importance, decreasing = TRUE),] ## this is to give an overview and choose cutoff
                  #picking up important with cutoff
                  BB1_sort_cutoff<-filter(BB1, Layer_importance >= 2) ## number is your importance cutoff

                        #picking up OTUs from full otutable
                        signOTU<-as.matrix(BB1_sort_cutoff$Row.names) #vector with OTUs
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = ".")
                        #my_taxa
                        tax$OTU <-rownames(tax)
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = "=")
                        signmy_taxa<-filter(tax, OTU %in% signOTU) 
                        signmy_taxa<-as.data.frame(signmy_taxa)
                        taxa_keep<-signmy_taxa$OTU
return(taxa_keep)
}

Run Baruta algorith with full data set. We are looking for taxa that explain differences among sites.

baruta_Site_indicators<-baruta_phyloseq(physeq = bb16S_bac_ML, GroupingVar = "Site")
baruta_Site_indicators_ps_object<-subset_taxa(bb16S_bac_ML, taxa_names(bb16S_bac_ML) %in% baruta_Site_indicators)

https://jokergoo.github.io/ComplexHeatmap-reference/book/ Heatmap to visualize differnces with all samples together.

Based off of initial figure, it looks like there is a strong grouping based on Baruta ID’ed taxa, and there is pattern for every other row. This indicates differences in time points within each sites. We might want to split by time point and then repeat analysis.

Since we observe what appears to be a time point signature, we can split by time point and rerun the Baruta algorithm.

#pre-treatment samples
bb16S_bac_ML_A<-subset_samples(bb16S_bac_ML, Timepoint == "A")
#run algo and subset phyloseq object
baruta_Site_indicators_A<-baruta_phyloseq(physeq = bb16S_bac_ML_A, GroupingVar = "Site")

After 18 iterations, +7.1 secs: 
 confirmed 3 attributes: OTU_21, OTU_5, OTU_8;
 rejected 2103 attributes: OTU_1000, OTU_1001, OTU_1002, OTU_1003, OTU_1004 and 2098 more;
 still have 151 attributes left.

After 23 iterations, +7.3 secs: 
 confirmed 9 attributes: OTU_1, OTU_117, OTU_197, OTU_2, OTU_20 and 4 more;
 rejected 16 attributes: OTU_10, OTU_105, OTU_109, OTU_119, OTU_137 and 11 more;
 still have 126 attributes left.

After 27 iterations, +7.4 secs: 
 confirmed 8 attributes: OTU_108, OTU_126, OTU_190, OTU_208, OTU_277 and 3 more;
 rejected 6 attributes: OTU_123, OTU_223, OTU_327, OTU_33, OTU_46 and 1 more;
 still have 112 attributes left.

After 30 iterations, +7.5 secs: 
 confirmed 1 attribute: OTU_80;
 rejected 9 attributes: OTU_154, OTU_186, OTU_19, OTU_22, OTU_236 and 4 more;
 still have 102 attributes left.

After 34 iterations, +7.7 secs: 
 confirmed 7 attributes: OTU_14, OTU_156, OTU_18, OTU_215, OTU_260 and 2 more;
 still have 95 attributes left.

After 37 iterations, +7.7 secs: 
 confirmed 2 attributes: OTU_196, OTU_451;
 rejected 1 attribute: OTU_130;
 still have 92 attributes left.

After 40 iterations, +7.8 secs: 
 confirmed 5 attributes: OTU_100, OTU_3, OTU_47, OTU_60, OTU_72;
 rejected 6 attributes: OTU_288, OTU_304, OTU_337, OTU_479, OTU_596 and 1 more;
 still have 81 attributes left.

After 44 iterations, +8 secs: 
 confirmed 4 attributes: OTU_127, OTU_143, OTU_16, OTU_24;
 rejected 1 attribute: OTU_163;
 still have 76 attributes left.

After 47 iterations, +8.1 secs: 
 confirmed 1 attribute: OTU_49;
 rejected 2 attributes: OTU_158, OTU_569;
 still have 73 attributes left.

After 50 iterations, +8.2 secs: 
 confirmed 5 attributes: OTU_110, OTU_144, OTU_193, OTU_62, OTU_83;
 rejected 1 attribute: OTU_220;
 still have 67 attributes left.

After 53 iterations, +8.2 secs: 
 confirmed 2 attributes: OTU_132, OTU_177;
 rejected 2 attributes: OTU_128, OTU_227;
 still have 63 attributes left.

After 55 iterations, +8.3 secs: 
 confirmed 1 attribute: OTU_85;
 still have 62 attributes left.

After 58 iterations, +8.4 secs: 
 rejected 1 attribute: OTU_180;
 still have 61 attributes left.

After 61 iterations, +8.5 secs: 
 confirmed 3 attributes: OTU_116, OTU_149, OTU_79;
 still have 58 attributes left.

After 64 iterations, +8.6 secs: 
 rejected 1 attribute: OTU_51;
 still have 57 attributes left.

After 67 iterations, +8.7 secs: 
 confirmed 1 attribute: OTU_15;
 rejected 1 attribute: OTU_125;
 still have 55 attributes left.

After 70 iterations, +8.8 secs: 
 confirmed 2 attributes: OTU_164, OTU_41;
 rejected 1 attribute: OTU_189;
 still have 52 attributes left.

After 72 iterations, +8.8 secs: 
 rejected 3 attributes: OTU_218, OTU_86, OTU_98;
 still have 49 attributes left.

After 88 iterations, +9.3 secs: 
 confirmed 1 attribute: OTU_7;
 still have 48 attributes left.

After 91 iterations, +9.4 secs: 
 confirmed 1 attribute: OTU_104;
 still have 47 attributes left.

After 101 iterations, +9.7 secs: 
 confirmed 1 attribute: OTU_35;
 still have 46 attributes left.

After 109 iterations, +9.9 secs: 
 rejected 1 attribute: OTU_150;
 still have 45 attributes left.

After 112 iterations, +10 secs: 
 confirmed 1 attribute: OTU_37;
 still have 44 attributes left.

After 119 iterations, +10 secs: 
 confirmed 1 attribute: OTU_231;
 still have 43 attributes left.

After 122 iterations, +10 secs: 
 confirmed 1 attribute: OTU_76;
 still have 42 attributes left.

After 132 iterations, +11 secs: 
 rejected 1 attribute: OTU_547;
 still have 41 attributes left.

After 176 iterations, +12 secs: 
 confirmed 1 attribute: OTU_17;
 still have 40 attributes left.

After 185 iterations, +12 secs: 
 confirmed 1 attribute: OTU_91;
 still have 39 attributes left.

After 195 iterations, +12 secs: 
 confirmed 1 attribute: OTU_211;
 still have 38 attributes left.

After 223 iterations, +13 secs: 
 rejected 1 attribute: OTU_358;
 still have 37 attributes left.

After 230 iterations, +13 secs: 
 rejected 1 attribute: OTU_293;
 still have 36 attributes left.

After 235 iterations, +13 secs: 
 rejected 1 attribute: OTU_78;
 still have 35 attributes left.

After 244 iterations, +14 secs: 
 confirmed 1 attribute: OTU_75;
 still have 34 attributes left.

After 256 iterations, +14 secs: 
 confirmed 1 attribute: OTU_65;
 still have 33 attributes left.

After 286 iterations, +15 secs: 
 confirmed 2 attributes: OTU_182, OTU_564;
 still have 31 attributes left.

After 291 iterations, +15 secs: 
 confirmed 1 attribute: OTU_66;
 still have 30 attributes left.

After 316 iterations, +16 secs: 
 rejected 1 attribute: OTU_204;
 still have 29 attributes left.

After 327 iterations, +16 secs: 
 rejected 1 attribute: OTU_289;
 still have 28 attributes left.

After 366 iterations, +17 secs: 
 rejected 1 attribute: OTU_102;
 still have 27 attributes left.

After 373 iterations, +17 secs: 
 rejected 1 attribute: OTU_292;
 still have 26 attributes left.

After 400 iterations, +18 secs: 
 confirmed 1 attribute: OTU_275;
 still have 25 attributes left.

After 404 iterations, +18 secs: 
 confirmed 1 attribute: OTU_298;
 still have 24 attributes left.

After 429 iterations, +18 secs: 
 confirmed 1 attribute: OTU_315;
 still have 23 attributes left.

After 440 iterations, +19 secs: 
 confirmed 1 attribute: OTU_270;
 still have 22 attributes left.

After 447 iterations, +19 secs: 
 confirmed 1 attribute: OTU_192;
 still have 21 attributes left.

After 467 iterations, +20 secs: 
 confirmed 1 attribute: OTU_448;
 still have 20 attributes left.

After 483 iterations, +20 secs: 
 rejected 1 attribute: OTU_361;
 still have 19 attributes left.

After 509 iterations, +21 secs: 
 confirmed 1 attribute: OTU_256;
 still have 18 attributes left.

After 521 iterations, +21 secs: 
 confirmed 1 attribute: OTU_226;
 still have 17 attributes left.

After 569 iterations, +22 secs: 
 confirmed 1 attribute: OTU_253;
 still have 16 attributes left.

After 626 iterations, +24 secs: 
 confirmed 1 attribute: OTU_162;
 still have 15 attributes left.

After 657 iterations, +24 secs: 
 confirmed 1 attribute: OTU_203;
 still have 14 attributes left.

After 720 iterations, +26 secs: 
 rejected 1 attribute: OTU_133;
 still have 13 attributes left.

After 914 iterations, +31 secs: 
 rejected 1 attribute: OTU_408;
 still have 12 attributes left.

After 939 iterations, +32 secs: 
 confirmed 1 attribute: OTU_131;
 still have 11 attributes left.

baruta_Site_indicators_ps_object_A<-subset_taxa(bb16S_bac_ML_A, taxa_names(bb16S_bac_ML_A) %in% baruta_Site_indicators_A)
# pull out taxonomic info
baruta_Site_indicators_ps_object_tax_A<-tax_table(baruta_Site_indicators_ps_object_A)
baruta_Site_indicators_ps_object_tax_A<-data.frame(baruta_Site_indicators_ps_object_tax_A)
#rename NA as unknown
baruta_Site_indicators_ps_object_tax_A$Class[is.na(baruta_Site_indicators_ps_object_tax_A$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_A$Family[is.na(baruta_Site_indicators_ps_object_tax_A$Family)]<-"Unknown"
#pull out taxonomy table for naming purposes
baruta_Site_indicators_ps_object_A@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_A)

Warning in .local(object) :
  Coercing from data.frame class to character matrix 
prior to building taxonomyTable. 
This could introduce artifacts. 
Check your taxonomyTable, or coerce to matrix manually.

fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_A)[,3]))

#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_A)[,4], row_dend_reorder =TRUE, col= col_fun)

#heatmap with otuIDs
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_A), row_dend_reorder =TRUE, col= col_fun)

Strong grouping of sites based on Baruta ID’ed taxa at the first sampling.

Again, strong groupings of sites by Baruta ID’ed taxa at time point 2.

Next, we will look for taxa that consistently explain the differences among treatments following application. Again we will use the Baruta algorithm, but this time, we will specify the grouping variable as Treatment.

Merge by treatments to see if there are consistencies across treatments. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment")
sample_data(bb16S_bac_ML_B_merged)$Treatment<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment")
#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below because no taxa are identified. 

baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)
#extract taxonomy and rename NA to unknown. 
baruta_Site_indicators_ps_object_tax_B_merged<-tax_table(baruta_Site_indicators_ps_object_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged<-data.frame(baruta_Site_indicators_ps_object_tax_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged$Class[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B_merged$Family[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B_merged@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B_merged)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B_merged)[,3]))
#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Merge by treatment site and rerun algo. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment_Site")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

sample_data(bb16S_bac_ML_B_merged)$Treatment_Site<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment_Site")

After 17 iterations, +3.2 secs: 
 rejected 1304 attributes: OTU_1, OTU_10, OTU_100, OTU_1002, OTU_1006 and 1299 more;
 no more attributes left.

#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below becasue no taxa are identified. 
baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)

Error in dimnames(x) <- dn : 
  length of 'dimnames' [1] not equal to array extent

Exploratory tree of indicators from initial Baruta and then split by timpeoints.

plot_tree(baruta_Site_indicators_ps_object, color="Site", shape = "Timepoint", label.tips = "taxa_names")

plot_tree(baruta_Site_indicators_ps_object_A, color="Site", shape = "Treatment", label.tips = "taxa_names")

plot_tree(baruta_Site_indicators_ps_object_B, color="Site", shape = "Treatment", label.tips = "taxa_names")

Shared and unique taxa among sites

unique_taxa(phyloseq_obj = bb16S_bac_rarefy, treatment = "Treatment")

$CF
  [1] "OTU_61"   "OTU_114"  "OTU_141"  "OTU_142"  "OTU_148"  "OTU_166"  "OTU_187" 
  [8] "OTU_237"  "OTU_240"  "OTU_255"  "OTU_311"  "OTU_326"  "OTU_384"  "OTU_390" 
 [15] "OTU_399"  "OTU_411"  "OTU_425"  "OTU_433"  "OTU_437"  "OTU_447"  "OTU_459" 
 [22] "OTU_469"  "OTU_477"  "OTU_487"  "OTU_492"  "OTU_500"  "OTU_503"  "OTU_513" 
 [29] "OTU_516"  "OTU_532"  "OTU_555"  "OTU_557"  "OTU_558"  "OTU_560"  "OTU_561" 
 [36] "OTU_565"  "OTU_566"  "OTU_583"  "OTU_586"  "OTU_591"  "OTU_601"  "OTU_609" 
 [43] "OTU_615"  "OTU_616"  "OTU_624"  "OTU_632"  "OTU_654"  "OTU_664"  "OTU_669" 
 [50] "OTU_678"  "OTU_679"  "OTU_690"  "OTU_696"  "OTU_697"  "OTU_704"  "OTU_705" 
 [57] "OTU_711"  "OTU_715"  "OTU_717"  "OTU_720"  "OTU_727"  "OTU_737"  "OTU_738" 
 [64] "OTU_744"  "OTU_748"  "OTU_750"  "OTU_755"  "OTU_761"  "OTU_762"  "OTU_765" 
 [71] "OTU_772"  "OTU_773"  "OTU_774"  "OTU_775"  "OTU_779"  "OTU_789"  "OTU_790" 
 [78] "OTU_791"  "OTU_803"  "OTU_809"  "OTU_817"  "OTU_818"  "OTU_820"  "OTU_827" 
 [85] "OTU_829"  "OTU_832"  "OTU_838"  "OTU_842"  "OTU_849"  "OTU_854"  "OTU_856" 
 [92] "OTU_863"  "OTU_868"  "OTU_869"  "OTU_872"  "OTU_873"  "OTU_874"  "OTU_876" 
 [99] "OTU_880"  "OTU_881"  "OTU_882"  "OTU_883"  "OTU_892"  "OTU_894"  "OTU_895" 
[106] "OTU_899"  "OTU_903"  "OTU_904"  "OTU_905"  "OTU_908"  "OTU_920"  "OTU_921" 
[113] "OTU_924"  "OTU_929"  "OTU_930"  "OTU_935"  "OTU_936"  "OTU_937"  "OTU_945" 
[120] "OTU_946"  "OTU_947"  "OTU_949"  "OTU_950"  "OTU_954"  "OTU_957"  "OTU_958" 
[127] "OTU_959"  "OTU_960"  "OTU_961"  "OTU_962"  "OTU_963"  "OTU_964"  "OTU_965" 
[134] "OTU_970"  "OTU_972"  "OTU_973"  "OTU_978"  "OTU_979"  "OTU_980"  "OTU_981" 
[141] "OTU_987"  "OTU_990"  "OTU_991"  "OTU_992"  "OTU_998"  "OTU_999"  "OTU_1000"
[148] "OTU_1005" "OTU_1011" "OTU_1016" "OTU_1019" "OTU_1020" "OTU_1029" "OTU_1041"
[155] "OTU_1045" "OTU_1046" "OTU_1054" "OTU_1055" "OTU_1059" "OTU_1060" "OTU_1061"
[162] "OTU_1070" "OTU_1071" "OTU_1072" "OTU_1073" "OTU_1080" "OTU_1083" "OTU_1085"
[169] "OTU_1086" "OTU_1090" "OTU_1091" "OTU_1092" "OTU_1093" "OTU_1098" "OTU_1102"
[176] "OTU_1103" "OTU_1106" "OTU_1107" "OTU_1108" "OTU_1110" "OTU_1112" "OTU_1113"
[183] "OTU_1118" "OTU_1119" "OTU_1120" "OTU_1132" "OTU_1136" "OTU_1137" "OTU_1138"
[190] "OTU_1142" "OTU_1143" "OTU_1148" "OTU_1152" "OTU_1153" "OTU_1154" "OTU_1155"
[197] "OTU_1156" "OTU_1169" "OTU_1170" "OTU_1171" "OTU_1172" "OTU_1173" "OTU_1179"
[204] "OTU_1180" "OTU_1181" "OTU_1185" "OTU_1186" "OTU_1187" "OTU_1190" "OTU_1191"
[211] "OTU_1192" "OTU_1198" "OTU_1199" "OTU_1209" "OTU_1210" "OTU_1211" "OTU_1212"
[218] "OTU_1213" "OTU_1215" "OTU_1216" "OTU_1220" "OTU_1221" "OTU_1226" "OTU_1227"
[225] "OTU_1228" "OTU_1229" "OTU_1236" "OTU_1237" "OTU_1251" "OTU_1256" "OTU_1257"
[232] "OTU_1259" "OTU_1263" "OTU_1265" "OTU_1266" "OTU_1268" "OTU_1269" "OTU_1270"
[239] "OTU_1273" "OTU_1276" "OTU_1286" "OTU_1293" "OTU_1294" "OTU_1295" "OTU_1301"
[246] "OTU_1302" "OTU_1304" "OTU_1306" "OTU_1307" "OTU_1313" "OTU_1314" "OTU_1315"
[253] "OTU_1323" "OTU_1327" "OTU_1328" "OTU_1337" "OTU_1338" "OTU_1340" "OTU_1341"
[260] "OTU_1349" "OTU_1360" "OTU_1362" "OTU_1363" "OTU_1371" "OTU_1372" "OTU_1373"
[267] "OTU_1374" "OTU_1375" "OTU_1395" "OTU_1396" "OTU_1397" "OTU_1399" "OTU_1409"
[274] "OTU_1415" "OTU_1416" "OTU_1417" "OTU_1418" "OTU_1422" "OTU_1423" "OTU_1424"
[281] "OTU_1425" "OTU_1426" "OTU_1433" "OTU_1434" "OTU_1435" "OTU_1436" "OTU_1437"
[288] "OTU_1451" "OTU_1452" "OTU_1453" "OTU_1454" "OTU_1461" "OTU_1462" "OTU_1463"
[295] "OTU_1464" "OTU_1465" "OTU_1468" "OTU_1469" "OTU_1470" "OTU_1471" "OTU_1478"
[302] "OTU_1484" "OTU_1486" "OTU_1487" "OTU_1488" "OTU_1489" "OTU_1490" "OTU_1491"
[309] "OTU_1500" "OTU_1501" "OTU_1502" "OTU_1503" "OTU_1516" "OTU_1529" "OTU_1530"
[316] "OTU_1531" "OTU_1532" "OTU_1535" "OTU_1542" "OTU_1543" "OTU_1544" "OTU_1550"
[323] "OTU_1551" "OTU_1552" "OTU_1553" "OTU_1565" "OTU_1566" "OTU_1567" "OTU_1579"
[330] "OTU_1580" "OTU_1581" "OTU_1582" "OTU_1583" "OTU_1584" "OTU_1585" "OTU_1586"
[337] "OTU_1605" "OTU_1606" "OTU_1607" "OTU_1630" "OTU_1632" "OTU_1633" "OTU_1643"
[344] "OTU_1644" "OTU_1645" "OTU_1647" "OTU_1648" "OTU_1649" "OTU_1659" "OTU_1660"
[351] "OTU_1661" "OTU_1662" "OTU_1663" "OTU_1672" "OTU_1673" "OTU_1675" "OTU_1676"
[358] "OTU_1683" "OTU_1686" "OTU_1688" "OTU_1689" "OTU_1690" "OTU_1720" "OTU_1721"
[365] "OTU_1722" "OTU_1723" "OTU_1748" "OTU_1749" "OTU_1750" "OTU_1751" "OTU_1752"
[372] "OTU_1753" "OTU_1754" "OTU_1755" "OTU_1757" "OTU_1771" "OTU_1772" "OTU_1773"
[379] "OTU_1784" "OTU_1785" "OTU_1786" "OTU_1788" "OTU_1789" "OTU_1801" "OTU_1803"
[386] "OTU_1805" "OTU_1807" "OTU_1809" "OTU_1822" "OTU_1824" "OTU_1825" "OTU_1826"
[393] "OTU_1827" "OTU_1829" "OTU_1831" "OTU_1833" "OTU_1837" "OTU_1838" "OTU_1839"
[400] "OTU_1892" "OTU_1894" "OTU_1895" "OTU_1928" "OTU_1929" "OTU_1930" "OTU_1931"
[407] "OTU_1932" "OTU_1933" "OTU_1934" "OTU_1936" "OTU_1937" "OTU_1938" "OTU_1941"
[414] "OTU_1952" "OTU_1953" "OTU_1968" "OTU_1969" "OTU_1970" "OTU_1971" "OTU_1972"
[421] "OTU_1973" "OTU_1974" "OTU_1975" "OTU_1976" "OTU_1977" "OTU_1978" "OTU_1981"
[428] "OTU_1982" "OTU_1985" "OTU_2006" "OTU_2007" "OTU_2008" "OTU_2009" "OTU_2010"
[435] "OTU_2011" "OTU_2012" "OTU_2014" "OTU_2015" "OTU_2016" "OTU_2017" "OTU_2018"
[442] "OTU_2019" "OTU_2038" "OTU_2041" "OTU_2042" "OTU_2045" "OTU_2046" "OTU_2048"
[449] "OTU_2050" "OTU_2051" "OTU_2054" "OTU_2055" "OTU_2118" "OTU_2121" "OTU_2122"
[456] "OTU_2124" "OTU_2126" "OTU_2161" "OTU_2163" "OTU_2164" "OTU_2168" "OTU_2170"
[463] "OTU_2193" "OTU_2197" "OTU_2198" "OTU_2199" "OTU_2201" "OTU_2202" "OTU_2203"
[470] "OTU_2204" "OTU_2205" "OTU_2227" "OTU_2231" "OTU_2232" "OTU_2233" "OTU_2234"
[477] "OTU_2235" "OTU_2236" "OTU_2238" "OTU_2240" "OTU_2241" "OTU_2242" "OTU_2243"
[484] "OTU_2245" "OTU_2246" "OTU_2279" "OTU_2280" "OTU_2283" "OTU_2287" "OTU_2289"
[491] "OTU_2291"

$CON
  [1] "OTU_73"   "OTU_94"   "OTU_212"  "OTU_252"  "OTU_344"  "OTU_352"  "OTU_369" 
  [8] "OTU_404"  "OTU_450"  "OTU_458"  "OTU_461"  "OTU_464"  "OTU_476"  "OTU_481" 
 [15] "OTU_489"  "OTU_510"  "OTU_525"  "OTU_526"  "OTU_534"  "OTU_537"  "OTU_542" 
 [22] "OTU_548"  "OTU_549"  "OTU_570"  "OTU_579"  "OTU_581"  "OTU_584"  "OTU_587" 
 [29] "OTU_593"  "OTU_610"  "OTU_617"  "OTU_621"  "OTU_627"  "OTU_638"  "OTU_644" 
 [36] "OTU_651"  "OTU_653"  "OTU_655"  "OTU_659"  "OTU_671"  "OTU_675"  "OTU_685" 
 [43] "OTU_686"  "OTU_691"  "OTU_707"  "OTU_713"  "OTU_718"  "OTU_719"  "OTU_721" 
 [50] "OTU_729"  "OTU_730"  "OTU_733"  "OTU_734"  "OTU_742"  "OTU_743"  "OTU_751" 
 [57] "OTU_757"  "OTU_763"  "OTU_766"  "OTU_776"  "OTU_785"  "OTU_787"  "OTU_794" 
 [64] "OTU_795"  "OTU_797"  "OTU_798"  "OTU_799"  "OTU_802"  "OTU_824"  "OTU_825" 
 [71] "OTU_830"  "OTU_831"  "OTU_835"  "OTU_839"  "OTU_840"  "OTU_844"  "OTU_846" 
 [78] "OTU_851"  "OTU_852"  "OTU_858"  "OTU_859"  "OTU_864"  "OTU_877"  "OTU_887" 
 [85] "OTU_888"  "OTU_890"  "OTU_898"  "OTU_900"  "OTU_915"  "OTU_919"  "OTU_922" 
 [92] "OTU_925"  "OTU_926"  "OTU_931"  "OTU_933"  "OTU_940"  "OTU_941"  "OTU_942" 
 [99] "OTU_943"  "OTU_948"  "OTU_952"  "OTU_956"  "OTU_966"  "OTU_967"  "OTU_974" 
[106] "OTU_975"  "OTU_976"  "OTU_983"  "OTU_984"  "OTU_985"  "OTU_988"  "OTU_993" 
[113] "OTU_997"  "OTU_1001" "OTU_1002" "OTU_1007" "OTU_1008" "OTU_1012" "OTU_1013"
[120] "OTU_1014" "OTU_1017" "OTU_1018" "OTU_1021" "OTU_1022" "OTU_1023" "OTU_1024"
[127] "OTU_1025" "OTU_1030" "OTU_1032" "OTU_1033" "OTU_1037" "OTU_1038" "OTU_1042"
[134] "OTU_1047" "OTU_1048" "OTU_1051" "OTU_1053" "OTU_1056" "OTU_1058" "OTU_1063"
[141] "OTU_1065" "OTU_1066" "OTU_1068" "OTU_1069" "OTU_1075" "OTU_1076" "OTU_1078"
[148] "OTU_1082" "OTU_1087" "OTU_1094" "OTU_1099" "OTU_1121" "OTU_1122" "OTU_1123"
[155] "OTU_1125" "OTU_1140" "OTU_1144" "OTU_1145" "OTU_1146" "OTU_1149" "OTU_1150"
[162] "OTU_1157" "OTU_1167" "OTU_1168" "OTU_1174" "OTU_1175" "OTU_1182" "OTU_1183"
[169] "OTU_1184" "OTU_1188" "OTU_1193" "OTU_1194" "OTU_1195" "OTU_1200" "OTU_1201"
[176] "OTU_1206" "OTU_1207" "OTU_1217" "OTU_1218" "OTU_1222" "OTU_1225" "OTU_1231"
[183] "OTU_1232" "OTU_1233" "OTU_1238" "OTU_1239" "OTU_1240" "OTU_1242" "OTU_1243"
[190] "OTU_1247" "OTU_1253" "OTU_1260" "OTU_1261" "OTU_1264" "OTU_1267" "OTU_1278"
[197] "OTU_1279" "OTU_1280" "OTU_1288" "OTU_1289" "OTU_1290" "OTU_1291" "OTU_1292"
[204] "OTU_1296" "OTU_1297" "OTU_1298" "OTU_1308" "OTU_1309" "OTU_1310" "OTU_1317"
[211] "OTU_1319" "OTU_1324" "OTU_1331" "OTU_1332" "OTU_1333" "OTU_1334" "OTU_1342"
[218] "OTU_1343" "OTU_1344" "OTU_1350" "OTU_1351" "OTU_1352" "OTU_1364" "OTU_1365"
[225] "OTU_1366" "OTU_1367" "OTU_1368" "OTU_1369" "OTU_1376" "OTU_1377" "OTU_1378"
[232] "OTU_1379" "OTU_1381" "OTU_1389" "OTU_1390" "OTU_1391" "OTU_1400" "OTU_1401"
[239] "OTU_1402" "OTU_1403" "OTU_1410" "OTU_1411" "OTU_1412" "OTU_1414" "OTU_1427"
[246] "OTU_1428" "OTU_1438" "OTU_1441" "OTU_1443" "OTU_1444" "OTU_1445" "OTU_1446"
[253] "OTU_1455" "OTU_1456" "OTU_1457" "OTU_1458" "OTU_1459" "OTU_1460" "OTU_1472"
[260] "OTU_1473" "OTU_1479" "OTU_1485" "OTU_1492" "OTU_1493" "OTU_1494" "OTU_1495"
[267] "OTU_1504" "OTU_1506" "OTU_1507" "OTU_1518" "OTU_1519" "OTU_1520" "OTU_1521"
[274] "OTU_1537" "OTU_1538" "OTU_1539" "OTU_1545" "OTU_1546" "OTU_1554" "OTU_1555"
[281] "OTU_1557" "OTU_1558" "OTU_1559" "OTU_1560" "OTU_1568" "OTU_1569" "OTU_1570"
[288] "OTU_1571" "OTU_1572" "OTU_1588" "OTU_1589" "OTU_1590" "OTU_1591" "OTU_1592"
[295] "OTU_1593" "OTU_1595" "OTU_1596" "OTU_1597" "OTU_1599" "OTU_1600" "OTU_1608"
[302] "OTU_1609" "OTU_1610" "OTU_1611" "OTU_1612" "OTU_1613" "OTU_1614" "OTU_1615"
[309] "OTU_1616" "OTU_1617" "OTU_1618" "OTU_1620" "OTU_1634" "OTU_1635" "OTU_1636"
[316] "OTU_1637" "OTU_1639" "OTU_1640" "OTU_1650" "OTU_1652" "OTU_1664" "OTU_1665"
[323] "OTU_1666" "OTU_1677" "OTU_1678" "OTU_1679" "OTU_1694" "OTU_1695" "OTU_1696"
[330] "OTU_1697" "OTU_1698" "OTU_1699" "OTU_1700" "OTU_1704" "OTU_1725" "OTU_1726"
[337] "OTU_1728" "OTU_1729" "OTU_1730" "OTU_1733" "OTU_1734" "OTU_1735" "OTU_1736"
[344] "OTU_1737" "OTU_1738" "OTU_1764" "OTU_1765" "OTU_1774" "OTU_1776" "OTU_1790"
[351] "OTU_1791" "OTU_1811" "OTU_1813" "OTU_1814" "OTU_1842" "OTU_1843" "OTU_1844"
[358] "OTU_1845" "OTU_1851" "OTU_1852" "OTU_1856" "OTU_1858" "OTU_1860" "OTU_1861"
[365] "OTU_1896" "OTU_1897" "OTU_1898" "OTU_1900" "OTU_1903" "OTU_1904" "OTU_1905"
[372] "OTU_1906" "OTU_1908" "OTU_1909" "OTU_1910" "OTU_1911" "OTU_1912" "OTU_1914"
[379] "OTU_1943" "OTU_1944" "OTU_1946" "OTU_1947" "OTU_1948" "OTU_1954" "OTU_1956"
[386] "OTU_1957" "OTU_1986" "OTU_1990" "OTU_1993" "OTU_1994" "OTU_2020" "OTU_2023"
[393] "OTU_2025" "OTU_2026" "OTU_2027" "OTU_2028" "OTU_2057" "OTU_2058" "OTU_2059"
[400] "OTU_2060" "OTU_2061" "OTU_2063" "OTU_2065" "OTU_2066" "OTU_2067" "OTU_2068"
[407] "OTU_2069" "OTU_2070" "OTU_2072" "OTU_2073" "OTU_2074" "OTU_2075" "OTU_2078"
[414] "OTU_2079" "OTU_2084" "OTU_2085" "OTU_2128" "OTU_2129" "OTU_2132" "OTU_2135"
[421] "OTU_2139" "OTU_2141" "OTU_2143" "OTU_2144" "OTU_2146" "OTU_2147" "OTU_2148"
[428] "OTU_2149" "OTU_2171" "OTU_2172" "OTU_2173" "OTU_2175" "OTU_2177" "OTU_2178"
[435] "OTU_2181" "OTU_2208" "OTU_2209" "OTU_2249" "OTU_2251" "OTU_2253" "OTU_2254"
[442] "OTU_2256" "OTU_2257" "OTU_2259" "OTU_2260" "OTU_2296" "OTU_2297" "OTU_2301"
[449] "OTU_2303" "OTU_2304" "OTU_2305" "OTU_2307" "OTU_2308" "OTU_2311"

$ORG
  [1] "OTU_48"   "OTU_232"  "OTU_244"  "OTU_318"  "OTU_371"  "OTU_376"  "OTU_414" 
  [8] "OTU_444"  "OTU_449"  "OTU_485"  "OTU_502"  "OTU_504"  "OTU_508"  "OTU_533" 
 [15] "OTU_539"  "OTU_540"  "OTU_544"  "OTU_588"  "OTU_595"  "OTU_607"  "OTU_623" 
 [22] "OTU_625"  "OTU_628"  "OTU_629"  "OTU_630"  "OTU_639"  "OTU_642"  "OTU_643" 
 [29] "OTU_647"  "OTU_650"  "OTU_663"  "OTU_695"  "OTU_703"  "OTU_709"  "OTU_731" 
 [36] "OTU_732"  "OTU_736"  "OTU_759"  "OTU_764"  "OTU_768"  "OTU_777"  "OTU_778" 
 [43] "OTU_782"  "OTU_792"  "OTU_801"  "OTU_805"  "OTU_806"  "OTU_807"  "OTU_811" 
 [50] "OTU_812"  "OTU_814"  "OTU_821"  "OTU_836"  "OTU_841"  "OTU_847"  "OTU_848" 
 [57] "OTU_853"  "OTU_855"  "OTU_857"  "OTU_867"  "OTU_870"  "OTU_875"  "OTU_886" 
 [64] "OTU_889"  "OTU_896"  "OTU_902"  "OTU_910"  "OTU_911"  "OTU_912"  "OTU_917" 
 [71] "OTU_927"  "OTU_928"  "OTU_938"  "OTU_939"  "OTU_969"  "OTU_986"  "OTU_995" 
 [78] "OTU_1003" "OTU_1010" "OTU_1015" "OTU_1026" "OTU_1027" "OTU_1028" "OTU_1034"
 [85] "OTU_1039" "OTU_1043" "OTU_1044" "OTU_1049" "OTU_1052" "OTU_1057" "OTU_1077"
 [92] "OTU_1084" "OTU_1088" "OTU_1095" "OTU_1096" "OTU_1097" "OTU_1101" "OTU_1104"
 [99] "OTU_1105" "OTU_1109" "OTU_1115" "OTU_1116" "OTU_1117" "OTU_1126" "OTU_1127"
[106] "OTU_1128" "OTU_1129" "OTU_1130" "OTU_1134" "OTU_1141" "OTU_1151" "OTU_1158"
[113] "OTU_1159" "OTU_1160" "OTU_1163" "OTU_1164" "OTU_1176" "OTU_1178" "OTU_1189"
[120] "OTU_1196" "OTU_1197" "OTU_1202" "OTU_1203" "OTU_1204" "OTU_1205" "OTU_1208"
[127] "OTU_1219" "OTU_1223" "OTU_1224" "OTU_1234" "OTU_1235" "OTU_1241" "OTU_1244"
[134] "OTU_1245" "OTU_1249" "OTU_1250" "OTU_1254" "OTU_1255" "OTU_1262" "OTU_1275"
[141] "OTU_1281" "OTU_1282" "OTU_1283" "OTU_1284" "OTU_1285" "OTU_1299" "OTU_1300"
[148] "OTU_1303" "OTU_1305" "OTU_1311" "OTU_1312" "OTU_1321" "OTU_1326" "OTU_1335"
[155] "OTU_1336" "OTU_1345" "OTU_1346" "OTU_1347" "OTU_1348" "OTU_1355" "OTU_1356"
[162] "OTU_1357" "OTU_1358" "OTU_1359" "OTU_1370" "OTU_1382" "OTU_1383" "OTU_1385"
[169] "OTU_1386" "OTU_1387" "OTU_1392" "OTU_1393" "OTU_1394" "OTU_1404" "OTU_1405"
[176] "OTU_1406" "OTU_1407" "OTU_1419" "OTU_1421" "OTU_1429" "OTU_1430" "OTU_1447"
[183] "OTU_1448" "OTU_1450" "OTU_1474" "OTU_1475" "OTU_1476" "OTU_1477" "OTU_1480"
[190] "OTU_1481" "OTU_1482" "OTU_1496" "OTU_1497" "OTU_1498" "OTU_1505" "OTU_1512"
[197] "OTU_1513" "OTU_1514" "OTU_1515" "OTU_1524" "OTU_1525" "OTU_1526" "OTU_1527"
[204] "OTU_1540" "OTU_1541" "OTU_1547" "OTU_1548" "OTU_1549" "OTU_1561" "OTU_1563"
[211] "OTU_1564" "OTU_1573" "OTU_1574" "OTU_1575" "OTU_1577" "OTU_1578" "OTU_1601"
[218] "OTU_1602" "OTU_1603" "OTU_1604" "OTU_1622" "OTU_1624" "OTU_1626" "OTU_1627"
[225] "OTU_1628" "OTU_1641" "OTU_1642" "OTU_1655" "OTU_1656" "OTU_1657" "OTU_1658"
[232] "OTU_1668" "OTU_1669" "OTU_1670" "OTU_1680" "OTU_1681" "OTU_1705" "OTU_1706"
[239] "OTU_1707" "OTU_1709" "OTU_1710" "OTU_1711" "OTU_1712" "OTU_1715" "OTU_1716"
[246] "OTU_1717" "OTU_1718" "OTU_1719" "OTU_1740" "OTU_1741" "OTU_1743" "OTU_1744"
[253] "OTU_1745" "OTU_1746" "OTU_1769" "OTU_1770" "OTU_1777" "OTU_1778" "OTU_1779"
[260] "OTU_1780" "OTU_1783" "OTU_1793" "OTU_1794" "OTU_1795" "OTU_1796" "OTU_1797"
[267] "OTU_1798" "OTU_1799" "OTU_1800" "OTU_1815" "OTU_1816" "OTU_1817" "OTU_1818"
[274] "OTU_1819" "OTU_1820" "OTU_1821" "OTU_1863" "OTU_1865" "OTU_1867" "OTU_1868"
[281] "OTU_1872" "OTU_1873" "OTU_1875" "OTU_1876" "OTU_1877" "OTU_1878" "OTU_1883"
[288] "OTU_1884" "OTU_1886" "OTU_1889" "OTU_1890" "OTU_1922" "OTU_1924" "OTU_1925"
[295] "OTU_1926" "OTU_1949" "OTU_1950" "OTU_1951" "OTU_1959" "OTU_1961" "OTU_1963"
[302] "OTU_1965" "OTU_1966" "OTU_1967" "OTU_1995" "OTU_1998" "OTU_1999" "OTU_2000"
[309] "OTU_2003" "OTU_2004" "OTU_2005" "OTU_2030" "OTU_2031" "OTU_2032" "OTU_2033"
[316] "OTU_2034" "OTU_2035" "OTU_2036" "OTU_2037" "OTU_2089" "OTU_2096" "OTU_2097"
[323] "OTU_2099" "OTU_2101" "OTU_2102" "OTU_2104" "OTU_2105" "OTU_2108" "OTU_2109"
[330] "OTU_2113" "OTU_2153" "OTU_2154" "OTU_2157" "OTU_2160" "OTU_2182" "OTU_2183"
[337] "OTU_2184" "OTU_2185" "OTU_2187" "OTU_2188" "OTU_2189" "OTU_2212" "OTU_2214"
[344] "OTU_2216" "OTU_2219" "OTU_2220" "OTU_2222" "OTU_2225" "OTU_2262" "OTU_2263"
[351] "OTU_2264" "OTU_2265" "OTU_2269" "OTU_2270" "OTU_2271" "OTU_2272" "OTU_2275"
[358] "OTU_2276" "OTU_2277" "OTU_2278" "OTU_2317" "OTU_2318" "OTU_2320" "OTU_2321"
[365] "OTU_2322" "OTU_2323"